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Statistical tools of uncertainty quantification can be used to assess the information content of 
measured observables with respect to present-day theoretical models; to estimate model errors and 
thereby improve predictive capability; to extrapolate beyond the regions reached by experiment; and 
to provide meaningful input to applications and planned measurements. To showcase new oppor¬ 
tunities offered by such tools, we make a rigorous analysis of theoretical statistical uncertainties in 
nuclear density functional theory using Bayesian inference methods. By considering the recent mass 
measurements from the Canadian Penning Trap at Argonne National Laboratory, we demonstrate 
how the Bayesian analysis and a direct least-squares optimization, combined with high-performance 
computing, can be used to assess the information content of the new data with respect to a model 
based on the Skyrme energy density functional approach. Employing the posterior probability distri¬ 
bution computed with a Gaussian process emulator, we apply the Bayesian framework to propagate 
theoretical statistical uncertainties in predictions of nuclear masses, two-neutron dripline, and fis¬ 
sion barriers. Overall, we find that the new mass measurements do not impose a constraint that is 
strong enough to lead to significant changes in the model parameters. The example discussed in this 
study sets the stage for quantifying and maximizing the impact of new measurements with respect 
to current modeling and guiding future experimental efforts, thus enhancing the experiment-theory 
cycle in the scientific method. 

PACS numbers: 21.10.Dr, 21.60.Jz, 24.75.+i, 02.30.Zz 


Introduction - Our understanding of heavy, complex 
nuclei lies at the heart of many basic science questions, 
such as chemical evolution, neutron star structure, syn¬ 
thesis of super heavy elements, mechanism of nuclear fis¬ 
sion, or search for the new Standard Model [T]; this 
knowledge is also crucial for societal applications [2]. In 
all those cases, reliable theoretical estimates of nuclear 
masses, low-lying excitations, electromagnetic strength, 
and nuclear reaction rates form essential inputs when di¬ 
rect experimental information is not available. 

For tackling complex nuclei theoretically, nuclear den¬ 
sity functional theory (DFT) is the microscopic tool of 
choice [3]. In recent years, largely because of algorith¬ 
mic developments and high-performance computing [4], 
DFT has taken great strides as a predictive theory that 
describes the properties of nuclei across the nuclear land¬ 
scape No consensus exists, however, on the form of 

the nuclear effective interaction or energy density func¬ 
tional (EDF), resulting in large systematic uncertainties. 
Moreover, nuclear EDFs are characterized by coupling 
constants that must be adjusted to experiment 0 Et 
HO]. The systematic calculation of uncertainties related 
to the determination of model parameters, as well as the 
propagation of these uncertainties in model prediction, 
has thus become a necessity 13 EHIS] (see also USD- 


Furthermore, as we enter the era of experiments with 
exotic nuclei at extremes of isospin, theory will play an 
increasingly important role in identifying scientific pri¬ 
orities of planned experimental campaigns. Conversely, 
as experiments extend current knowledge by providing 
information about the uncharted regions of the nuclear 
landscape, new methodologies become critical for evalu¬ 
ating the impact of these measurements on theory. 

From the viewpoint of statistics, determining the pa¬ 
rameters of a model given a set of experimental data mea¬ 
surements is an inverse problem m- Bayesian inference 
methods m are one of the most popular and powerful 
statistical approaches to inverse problems, with diverse 
applications in physics [T9j[20] (for recent nuclear physics 
applications, see, e.g.. Refs. |2lH2l])- In the Bayesian set¬ 
ting, model parameters are treated as random variables, 
and their uncertainty is characterized by their joint prob¬ 
ability distribution. Various techniques, often based on 
Monte Carlo simulations, have been developed to recon¬ 
struct this probability distribution from model prediction 
of experimental data. 

Objectives - In this work, we present the advanced ap¬ 
plication of Bayesian inference to global nuclear prop¬ 
erties using nuclear DFT. In particular, we use the 
Bayesian framework to quantify and propagate DFT sta- 
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tistical model uncertainties and to assess the information 
content of new data with respect to model developments. 
To this end, we study the impact of the recently reported 
mass measurements from the Canadian Penning Trap 
(CPT) mass spectrometer at Argonne National Labora¬ 
tory [28H3Q] on the Bayesian posterior probability distri¬ 
bution as well as the direct determination of EDF pa¬ 
rameters. The CPT dataset is unique in that it probes 
neutron-rich nuclei around hence, it can help im¬ 

prove our knowledge of isovector EDF properties and re¬ 
duce extrapolation uncertainties into the region of the 
astrophysical r-process. From the resulting posterior dis¬ 
tribution, we assess model uncertainties on observables, 
including the position of the two-neutron dripline and 
fission barrier heights of actinide nuclei. 

Method - Our theoretical framework is nuclear density 
functional theory with Skyrme EDFs. Pairing is mod¬ 
eled with a density-dependent pairing force and treated 
at the Hartree-Fock-Bogoliubov (HFB) level by using an 
approximate particle number projection with the Lipkin- 
Nogami method. We choose the unedfI parameteriza¬ 
tion of the Skyrme functional as our reference model [31] . 
This EDF is characterized by twelve parameters that 
were optimized on a set of binding energies for spherical 
and deformed nuclei, charge radii, odd-even mass differ¬ 
ences, and excitation energies of selected fission isomers 
(see Refs. [3TI433] for details of the model and the unedf 
EDF family). 

The quality of the functional is measured by a com¬ 
posite function. 
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where x denotes the set of model parameters, = 12 
is the number of model parameters, ut the number of 
different data types used in the fit {tit = 4 in our case), 
fit is the number of data points used for each data type, 
is the total number of data points, and dtj 
and ytj{x) are the experimental value and correspond¬ 
ing model prediction, respectively, for the jth data point 
of type t. For the unedf 1 functional, where rid = 115, 
computing the requires about 5 minutes of CPU time 
with over 800 cores in a multithreaded implementation of 
the DFT solver HFBTHO [34] . Monte Carlo simulations 
used to construct the posterior distribution may typically 
involve tens of thousands of such evaluations; even 
with current supercomputers, this cost is too high. We 
thus replace the DFT model ytj{x) with a Gaussian pro¬ 
cess (GP) response surface, allowing Monte Carlo-based 
Bayesian computation. 

The GP response surface is estimated within the en¬ 
compassing Bayesian formulation [35] by using an ensem¬ 
ble of DFT runs for each of the rid experimental nuclei 
used in 0- The ensemble is defined by a 200 x rix matrix 
of input settings distributed according to a space-filling 


Latin hypercube sample [36] over an -dimensional hy¬ 
perrectangle centered on the unedf 1 values. For each 
parameter, widths are determined according to the stan¬ 
dard deviations reported in Ref. [31], which were ob¬ 
tained through a covariance analysis that assumed a lin¬ 
ear approximation. The GP is controlled by a scaling 
parameter, as well as correlation parameters regulating 
the smoothness of the response surface in each of the rix 
parameter directions. 

The full posterior density includes a likelihood term 
for the experimental data based on Eq. 0 and the en¬ 
semble of training runs for the GP, the uniform prior for 
the model parameters x, and priors for the parameters 
that control the GP-based response surface; see Ref. m\ 
for a detailed description of the posterior density. We 
construct dependent samples from this distribution us¬ 
ing Markov chain Monte Carlo as detailed in [35], from 
which summaries such as 90% probability intervals and 
posterior means can be constructed. 

Results - Through Bayesian model calibration, we first 
obtained the posterior probability distribution for the 
UNEDF 1 parameter set, which provides a sense of how 
the set of fit observables of unedf 1 constrains the pa¬ 
rameters. In Fig. we show the univariate and bivari¬ 
ate marginal estimates of the posterior distribution. The 
blue-outlined regions give the 95% posterior probability 
region for the original UNEDF 1 parameters. We notice 
that the Bayesian approach is in agreement with esti¬ 
mates of uncertainties based on covariance analysis re¬ 
ported in Ref. m- In particular, most distributions are 
centered on the UNEDF 1 values, and the standard devia¬ 
tions extracted from the distribution are consistent with 
the 95% probability intervals. 

In a second step, we used our Bayesian formulation to 
evaluate the information content of the new mass mea¬ 
surements [28U3Q] . To this end, we modified the of 
Eq. 0 to include 17 new masses of neutron-rich even- 
even nuclei measured at the CPT; the experimental val¬ 
ues are listed in the supplemental material [38]. The 
GP response surface was again produced by using an 
augmented ensemble of 200 {rid + 17) DFT model eval¬ 
uations. The green-outlined regions in Fig. represent 
the same 95% posterior probability regions obtained with 
the inclusion of the Argonne mass measurements. With 
the exception of a few ill-constrained parameters (e.g., 
nuclear incompressibility and isovector surface coupling 
constant), the shift in the posterior is small for each pa¬ 
rameter. This suggests a weak impact of the additional 
data on our model. 

For comparison, we performed a direct reoptimization, 
independent of the GP response surface, of the unedf 1 
functional that includes the new CPT masses [39] . We re¬ 
fer to the reoptimized EDF parameter set as unedf Icpt; 
see supplemental material for parameter values [38] . The 
two parameterizations are similar. The largest relative 
difference, weighted by the standard deviations reported 
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FIG. 1. (Color online) Univariate and bivariate marginal es¬ 
timates of the posterior distribution for the 12-dimensional 
DFT parameter vector of the unedfI parameterization. The 
blue lines enclose an estimated 95% region for the poste¬ 
rior distribution found when only the original unedfI data 
are accounted for; the green-outlined regions represent the 
same region for the posterior distribution found when the 
new CPT mass measurements are included. The ranges 
of parameter variations are 0.155 < Pc < 0.165 (fm“^); 
-16.0 < E^^/A < -15.5 (MeV); 200 < < 240 

(MeV); 28 < < 30 (MeV); 20 < < 60 (MeV); 

0.8 < 1/M; < 1.2; -60 < < -40 (MeVfm^); -200 < 

^ (MeVfm^); -200 < UJ' < -150 (MeVfm^); 
-220 < < -180 (MeVfm^); -80 < < -60 

(MeVfm^); and -80 < < 0 (MeVfm^). 


TABLE I. Root-mean-square deviations for each of the types 
of data included in the unedf optimization. Masses and en¬ 
ergies are in MeV, radii in fm. 


Class 

unedfI 

unedfIcpt 

masses (def) 

0.721 

0.578 

masses (sph) 

1.461 

1.545 

radii 

0.022 

0.022 

odd-even staggering (n) 

0.023 

0.024 

odd-even staggering (p) 

0.079 

0.081 

fission isomer energies 

0.190 

0.316 

masses (CPT) 

1.064 

0.479 


in Ref. [31], is 0.6cr for the isovector surface coupling con¬ 
stant and surface symmetry energy. These quanti¬ 

ties have been difficult to constrain with the dataset used 
in the unedf protocol. Of interest, then, is the fact that 
the UNEDF 1 CPT value of is close to that of unedf2, 

which was also optimized to effective single-particle en¬ 
ergies [33] , known to be sensitive probes of surface prop¬ 
erties. Since the new dataset including CPT masses is 
more skewed toward neutron-rich nuclei, it may supply 


additional information about the shell structure above 
doubly magic ^^^Sn through a better determination of 
isovector coupling constants. 

Table [T] displays the root-mean-square deviation be¬ 
tween calculated and measured values for each type of 
data included in the optimization. We note that the in¬ 
clusion of the CPT mass measurements shifts the opti¬ 
mization priority, so that the new masses and deformed 
masses are reproduced more accurately, while predic¬ 
tions for fission isomers and spherical masses deteriorate 
slightly. The results in Table |l] are indicative of a small, 
additional constraint on the isovector coupling constants 
in unedfIcpx- 
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FIG. 2. (Golor online) Estimated theoretical error bars for 
the masses of the even-even nuclei measured in Refs. [28H30] . 
using the posterior for unedfI. Dark blue bands represent 
the 90% confidence bands obtained from the posterior; larger, 
light blue bands also account for model error; black bars show 
mass residuals. 

Equipped with the posterior distribution for the EDF 
parameters, we now turn to the propagation of statis¬ 
tical uncertainties for model predictions. We take the 
posterior distribution for the EDF parameters obtained 
by conditioning only on the UNEDF 1 measurements and 
propagate the distribution through the augmented GP- 
based emulator, producing prediction intervals for the 
new CPT mass measurements. These estimates are gen¬ 
uine holdout predictions since the new mass data were not 
used in determining the posterior distribution. Figure 
shows 90% prediction intervals (centered on the mean 
mass value of unedfI) for the new CPT masses. The 
dark blue band is the 90% interval for the uncertainty in 
the EDF model parameters; the light blue band also in¬ 
cludes uncertainty due to model error. The model error 
uncertainty was estimated from the difference between 
the posterior mean estimate and the actual mass mea¬ 
surements in the unedfI dataset. Separate estimates 
were made for spherical and deformed nuclei. These es¬ 
timated model-error standard deviations were assumed 
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to be appropriate for the new CPT mass measurements, 
producing this additional uncertainty. We observe that 
the experimentally measured values (black bars in Fig.[^ 
are generally within the 90% prediction interval. The es¬ 
timated uncertainty for the calculated masses is approxi¬ 
mately ±2 MeV, and slightly larger for the four spherical 
nuclei (the first four nuclei in the figure). This uncer¬ 
tainty is relatively large and in excellent agreement with 
the r.m.s. deviation for masses of even-even nuclei across 
the entire nuclear landscape, which is 1.9 MeV for UN- 
EDFl. 



FIG. 3. (Color online) Comparison between the two-neutron 
dripline predictions made with unedfI (solid line) and those 
made with unedfIcpt (dashed line). The 90% probability 
spread about the unedfI predictions is shown in grey. 

We now evaluate how the calculated model uncertain¬ 
ties impact predictions for important physical observ¬ 
ables. We first look at the position of the two-neutron 
dripline, which is especially important for our under¬ 
standing of nucleosynthesis in the r-process [40]. For a 
given element characterized by its proton number Z, the 
two-neutron dripline is defined as the point where the 
two-neutron separation energy becomes negative. We 
have performed an ensemble of calculations of nuclear 
binding energies for all even-even neutron-rich elements 
with 20 < Z < 100 over the Latin hypercube sample 
design of EDF parameter inputs, allowing yet another 
GP-based emulator to be constructed for these binding 
energies. 

Once the emulator is constructed, we propagate the 
posterior distribution of the model parameters (condi¬ 
tioning on either the unedfI or unedfIcpt datasets), 
producing uncertainty in the estimated dripline. With 
this Monte Carlo sample, we can estimate the posterior 
mode and 90% interval for the dripline for each value of 
Z. We explored the axial quadrupole potential energy 
surface of each nucleus to allow for deformed solutions. 
The results are presented in Fig. We observe that 


the inclusion of 17 new masses of neutron-rich nuclei in 
the optimization protocol did not impact the position of 
the dripline, since results with unedfI and unedfIcpt 
are practically indistinguishable. The predicted dripline 
is consistent with the results of large-scale DFT surveys 
013. Apart from the few closed-shell, waiting-point nu¬ 
clei, the uncertainty on the position of the dripline is 
on the order of 15 to 20 nucleons. This is comparable 
to statistical and systematic uncertainties obtained by 
comparing predictions made with different Skyrme func¬ 
tionals |6|. 



FIG. 4. (Golor online) Gomparison between the fission bar¬ 
rier predictions for made with unedfI (solid line) and 

those made with unedfIcpt (dashed line), together with the 
90% confidence interval (shaded grey area). The potential 
energy surface was obtained by following the lowest-energy 
static fission pathway in a four-dimensional collective space 
of axial and triaxial quadrupole, axial octupole, and axial 
hexadecapole mass moments. 

Another important application area of nuclear DFT 
is fission theory. In Fig. we show the potential en¬ 
ergy curve of ^^^Pu. This nucleus is representative of the 
actinide region and is often used as a theoretical bench¬ 
mark. Again, the results of unedfI and unedfIcpt are 
close. The large theoretical uncertainty in the predicted 
static fission barrier is worth noting; similar results were 
obtained in Ref. mi in the context of fission properties 
for r-process nuclei. Since a 1 MeV shift in the fission bar¬ 
rier translates into many orders of magnitude difference 
in the spontaneous fission half life, such results highlight 
the urgent need for better constraining the deformation 
properties of current EDFs. 

Conclusions - We have presented a comprehensive ap¬ 
plication of Bayesian inference techniques to the calcu¬ 
lation and propagation of theoretical statistical uncer¬ 
tainties in nuclear density functional theory. By using 
the recent, unique dataset of mass measurements from 
the CPT at Argonne National Laboratory, we showcase 
how the statistical tools of uncertainty quantification and 
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high-performance computing can be used to assess the 
information content of new data with respect to current 
models. Such analyses will become increasingly relevant 
for enhancing the feedback in the “observation-theory- 
prediction-experiment”- cycle of the scientific method at 
the eve of next-generation radioactive ion beam facilities 
and exascale computing. 

In the particular case studied in this work, we found 
that the impact of the new neutron-rich nuclei mass 
data on our DFT model is minor. The coupling con¬ 
stants of the earlier functional unedfI and of the new 
functional unedfIcpt, informed by the new data, are 
fairly close; hence, their predictions for the two-neutron 
dripline and fission barrier in ^^^Pu are practically iden¬ 
tical. Although the major theoretical statistical uncer¬ 
tainty in developments of the nuclear EDF comes from 
the poorly constrained isovector terms and the new data 
on neutron-rich nuclei are generally expected to reduce 
this uncertainty, the lack of a significant constraint from 
the new masses suggests that both the amount of new 
neutron-rich isotope data and the range of neutron excess 
probed, are not sufficiently large to impact our model ap¬ 
preciably. Moreover, because of their poor precision with 
respect to the existing data (see Table |I| , even the cur¬ 
rent, best-calibrated EDFs are not sensitive and flexible 
enough to fully take advantage of the new experimental 
information. 

By propagating theoretical errors, we found large 
model uncertainties in the predictions of the two-neutron 
dripline and the fission barrier in ^^^Pu. In this respect, 
we concur with the conclusions of Ref. [30] that existing 
mass models are insufficient for accurate r-process simu¬ 
lations. Clearly, accurate measurements for nuclei with 
even larger neutron excess, closer to the r-process path, 
are still needed in order to better inform theory. 

We note that the uncertainties discussed in this work 
are estimated statistically, reflecting parameter uncer¬ 
tainty and model misfit. The misfit error is most likely 
due to our lack of knowledge of the form of the nuclear 
EDF itself, and additional measurements will never re¬ 
duce this source of uncertainty. Adding physics that is 
missing in the current implementations of nuclear DFT 
is a major challenge for the field. A distinct and comple¬ 
mentary challenge is to develop tools that deliver uncer¬ 
tainty quantification for theoretical studies as well as for 
the assessment of new experimental data. The present 
work represents a step in this direction. 
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Experimental datasets 


UNEDF coupling constants 


The experimental dataset of unedfI contains Ud = 
115 data points [1, 2], which can be broken down into 
riT = 4 data types: rii = 75 nuclear masses (28 spherical 
and 47 deformed), n 2 = 28 r.m.s. proton radii, ns = S 
odd-even mass staggering differences (4 for neutrons and 
4 for protons), and 77,4 = 4 excitation energies of fis¬ 
sion isomers. Compared with unedfI, the unedfIcpt 
dataset contains 17 new masses of neutron-rich even-even 
nuclei measured by using the Canadian Penning Trap 
mass spectrometer and CARIBU facility [3-5] at Argonne 
National Laboratory. These new data are listed in Ta¬ 
ble I. 


TABLE 1. Experimental binding energies (rounded to the 
nearest 0.1 MeV) of the 17 even-even nuclei measured in 
Refs. [3-5] included in our analysis. 

Nucleus 

B (MeV) 

Ref. 


-1090.2 

[5] 


-1102.7 

[5] 


-1108.8 

[5] 

i34Te 

-1123.3 

[4] 

i36Te 

-1131.3 

[4] 

i38Te 

-1138.7 

[5] 

i40Te 

-1145.7 

[5] 


-1151.4 

[4] 


-1160.6 

[4] 


-1180.0 

[3] 


-1190.1 

[3] 

i«Ba 

-1199.4 

[3] 


-1208.5 

[3] 


-1219.4 

[3] 


-1230.0 

[3] 


-1291.8 

[4] 


-1302.9 

[4] 


Table II lists the coupling constants of the unedfO [1], 
unedfI [ 2 ], and unedfIcpt (this work) energy density 
functionals. 

TABLE 11. Coupling constants of the unedfO, unedfI, 
and unedfIcpt energy density functionals, pc is in fm“^; 

^-nd Lf^ are in MeV; 1/M* is dimen¬ 
sionless; and C^^'^ are in MeVfm^; and and Vq are 

in MeV fm^. 


Name 

unedfO 

unedfI 

unedfIcpt 

Pc 

0.1605 

0.1587 

0.1589 

E^^/A 

-16.0559 

-15.8000 

-15.8000 


230.0000 

220.0000 

220.0000 

NM 

'-t'sym 

30.5429 

28.9362 

29.3449 

r NM 
-^sym 

45.0804 

40.0149 

40.7144 

iim: 

0.9000 

0.9924 

0.9686 

QpAp 

-55.2606 

-45.1289 

-43.9801 

QpAp 

-55.6226 

-145.3178 

-114.2915 

d 

-170.3740 

-186.0655 

-182.2372 

U 

-199.2020 

-206.5796 

-203.9807 

J 

-79.5308 

-74.0264 

-72.4172 

J 

45.6302 

-35.6584 

-32.9206 
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